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ABSTRACT 

One-dimensional  stress  and  temperature  fields  in  a suddenly 
loaded  and/or  heated  semi-infinite  rod  of  nonlinear  thermoviscoelastic 
material  are  studied  using  coupled  thermomechanical  theory.  The 
transport  of  heat  is  governed  by  the  modified  Fourier  heat  conduction 
law,  and  thus  proceeds  by  wave  propagation  rather  than  by  diffusion. 

The  application  of  thermal  and  mechanical  disturbances  at  the  end  of 
the  rod  gives  rise  to  two  wave  fronts  along  which  these  disturbances 
propagate.  Field  solutions  for  the  stress  and  temperature  are 
obtained  by  numerical  integration  along  the  five  characteristics  of 
the  governing  equations,  and  results  are  presented  for  several  linear 
and  nonlinear  viscoelastic  models.  A closed  form  solution  is  also 
obtained  via  the  Laplace  transform  for  the  special  case  of  an  uncoupled 
thermal  wave. 
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specific  heat  at  constant  stress 
Young's  modulus 

integer  as  an  indication  of  position 

integer  for  the  condensation  of  computer  storage. 

Equation  (42) 

Isotropic  thermal  conductivity 
nondimens ional  material  constants 
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number  of  nonlinear  memory  integrals,  Equation  (2) 
number  of  components  of  strain,  Equation  (7) 
steady  creep  power.  Equation  (2) 
one-dimensional  heat  flux,  Equation  (1) 
transient  creep  powers,  Equation  (2) 
temperature 

constant  reference  temperature 
time 

n ’'"time  scale  for  nondimensionalization 
coupled  wave  speeds.  Equation  (6a) 
uncoupled  elastic  mechanical  wave  speed 
uncoupled  thermal  wave  speed 
particle  velocity 
space  coordinate 

coefficient  of  thermal  expansion 

positive  nondimensional  quantities  governing  the  wave 
speeds.  Equations  (5) 


i 
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one-dimensional  strain 


vi 


linear  elastic  strain,  Equation  (7b) 
steady  creep  strain.  Equation  (7c) 

e^.i-3 N transient  creep  strains.  Equation  (7d) 

et  thermal  strain,  Equation  (7e) 

6<*T-T  temperature  increment  relative  to  constant  reference 

temperature  Tq 

0q  input  temperature  discontinuity 

material  constants,  Equation  (2) 

p mass  density 

a one-dimensional  stress 

o input  stress  discontinuity 

o 

t relaxation  time  of  heat  conduction,  Equation  (1) 


Uj.J-1.2 

(") 


retardation  time  of  transient  creep,  Equation  (2) 
indicates  a discontinuity  across  the  leading  and  lagging 
wave  fronts  respectively 

indicates  nondimensional  variables  and  parameters 
(dropped  after  Equation  (9)) 

indicates  nondimensional  variables  for  the  uncoupled 


thermal  wave  problem 


The  classical  Fourier  heat  conduction  law  is  based  on  the  hypo- 
thesis that  the  heat  flux  is  linearly  proportional  to  the  temperature 
gradient,  and  predicts  an  infinite  thermal  wave  speed.  This  physi- 
cally anomalous  behavior  has  prompted  a number  of  proposed  modifications 
to  the  heat  conduction  law  which  can  predict  two  thermomechanically 
coupled  wave  speeds  (i.e.,  second  sound.)  The  second  sound  effect  has 
been  experimentally  observed  and  extensively  studied  in  solids  at  low 
temperatures  (e.g.,  see  Rogers  [1]),  and  it  has  been  noted  that 
second  sound  disturbances  experience  damping  effects.  From  the  point 
of  view  of  continuum  mechanics,  damping  may  result  from  an  additional 
term  in  a modified  heat  conduction  law,  a thermomechanical  coupling 
term  in  the  energy  equation,  and  from  viscoelastic  terms  in  the  consti- 
tutive relation.  The  relative  importance  of  these  three  damping 
mechanisms  may  vary  with  temperature,  but  they  are  in  general  all 
present.  It  is  our  primary  goal  to  demonstrate  that  even  if  all  these 
damping  mechanisms  are  included  in  a general  way  resulting  in  a very 
difficult  nonlinear  boundary  value  problem,  the  problem  may  still 
be  successfully  treated  numerically  through  the  use  of  the  method  of 
characteristics. 

Chester  [2]  employed  a modified  Fourier  heat  conduction  law  in 
one-dimension  as 
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12 

dt 


+ Q 


(1) 


TM’*' 


I 
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and  showed  that  the  temperature  obeys  a dissipative  wave  equation  in 
space  x and  time  t . In  (1)  Q is  the  heat  flux,  0 is  the 


2 


temperature  difference  T - Tq  relative  to  the  constant  reference 


temperature  T 


k is  the  thermal  conductivity  and  t i6  the  relax 


at ion  time.  Lord  and  Shulman  [3],  and  Achenbach  [4]  have  employed 
(1)  to  study  waves  in  linear  thermoelastic  materials,  and  recently 
Chang  and  Cozzarelli  [5]  have  extended  these  results  to  nonlinear 
thermoviscoelastlc  materials. 

Viscoelastic  constitutive  relations  may  be  developed  from  the 
laws  of  thermodynamics  by  the  energy  functional  approach  or  by  the 
state  variable  approach.  Both  approaches  yield  similar  results  as  dls 
cussed  by  Cost  [6].  In  [5]  a three-dimensional  single  integral  con- 
stitutive equation,  of  the  generalized  Kelvin  type,  was  obtained  via 
the  energy  functional  approach  for  nonlinear  thermoviscoelastlc 
materials.  For  one-dimensional  stress  o , strain  e and  temperature 
this  relation  may  be  written  as 


-Lt-t')  3|~|qi 


f 3|r!  r M ■ ' 'u. 1 

g+j(t-t')  sgn(o)dt ' +j  l (1-e  Ti  — sgn(o)dt ' 


+a6 

(2) 


where  sgn(o)  is  the  signum  function.  In  the  above  n is  the  steady 
creep  power  and  X is  the  steady  creep  parameter;  M is  the  number 
of  transient  creep  memory  integrals,  q^  are  the  transient  creep 
powers,  are  transient  creep  parameters  and  are  retardation 

times;  and  a is  the  coefficient  of  thermal  expansion.  This  kind  of 
single  integral  representation  is  quite  general  and  is  convenient  for 
the  development  of  analytical  and  numerical  solutions. 

A thermomechanically  coupled  linearized  energy  equation  identi- 
cal with  the  result  of  linear  thermoelasticity  was  also  obtained  in 
[5]  as 
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+ aT  I2-  + pC  ~ 
o 3t  a 3t 
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(3) 


where  p is  the  mass  density  and  is  the  specific  heat  at  constant 

stress.  In  obtaining  Equation  (3)  the  energy  dissipation  and  the  ther- 
mal kinetic  energy  (see  Kaliski  [7]  for  details)  have  been  dropped  as 
higher  order  terms. 

Section  2 contains  the  governing  equations  of  one-dimensional 
coupled  nonlinear  thermoviscoelastic  wave  propagation  with  finite 
thermal  wave  speed,  for  the  problem  of  a semi-infinite  rod  subjected 
to  sudden  thermal  and  mechanical  disturbance  at  its  end.  These  equa- 
tions include  the  previously  discussed  equations  (1-3)  plus  the  strain- 
displacement  relation  and  equation  of  motion  for  small  deformation 
theory.  Using  the  prescribed  initial  values  and  appropriate  reference 
quantities,  a nondimensional  form  of  the  governing  equations  is  also 
obtained.  In  section  3 the  method  of  characteristics  is  presented  in 
full  detail  for  the  problem  posed.  A closed  form  solution  based  on 
the  Laplace  transform  as  well  as  a numerical  solution  based  on  the 
method  of  characteristics  are  presented  in  section  4 for  the  special 
case  of  uncoupled  thermal  wave  propagation.  The  details  of  a numerical 
procedure  for  the  general  problem  employing  numerical  integration 
along  five  characteristics  in  a finite  difference  mesh  are  given  in 
section  5,  along  with  some  discussion  of  numerical  results  for  several 


linear  and  nonlinear  viscoelastic  models. 


The  final  section  contains 


a brief  summary. 


2.  GOVERNING  EQUATIONS  FOR  THE  SEMI-INFINITE  ROD 


Consider  a straight,  semi-infinite  rod  of  nonlinear  thermo- 
viscoelastic  material,  which  is  assumed  to  be  initially  stress  free 
and  at  uniform  temperature  Tq  . A constant  stress  oq  and  a con- 
stant temperature  increment  0q  are  suddenly  applied  at  the  end  of 
the  rod  (x  ■ 0) , and  two  thermomechanlcally  coupled  wave  fronts  are 
generated  along  the  positive  x axis.  We  designate  the  wave  speed 
of  the  leading  wave  by  and  that  of  the  lagging  wave  by  , 

where  and  V2  are  the  two  positive  roots  of  the  biquadratic 

equation  (see  [5]) 


<^)  - <1  + « + Y)(is^)2  + Y - 0 


(4) 


Here,  Vg  ■ /E/p  is  the  uncoupled  (a  m 0)  elastic  mechanical  wave 
speed,  and 


> 0 


Et(oC  -T  Ea  ) 
K o o 


T Ea 

6 • o 

pC  -T  Ea 
a o 
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(5b) 


are  non-negative  nondlmensional  material  constants,  where  6 ■ 0 
is  the  uncoupled  case.  Equation  (4)  then  yields 
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In  the  above  N is  the  number  of  components  of  strain,  repre- 
sents the  linear  elastic  strain,  Is  the  steady  creep  strain, 

for  1 - 3,...,N  represents  the  components  of  transient  creep 


strain,  and  et  Is  the  thermal  strain. 


For  Infinitesimal  deformation,  the  one  dimensional  strain- 
displacement  relation  and  the  equation  of  motion  are  given  by 
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where  v is  the  particle  velocity  in  the  x direction.  Equations 
(1,3, 7, 8),  with  the  appropriate  initial  and  boundary  conditions,  are 
the  complete  set  of  governing  equations. 


2(b).  Governing  Equations  in  Nondimensional  Form 
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from  which  it  also  follows  that 
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(9a) 


(9b) 


The  bar  as  used  in  Equations (9)  will  be  dropped  hereafter  for 
convenience  of  notation,  and  all  subsequent  variables  and  parameters 
will  be  nondimensional.  The  nondimensional  governing  equations  are 
then  written  as 


where  • Equations  (10)  are  N + 5 first  order  partial 

differential  equations  in  the  variables  e.(i“l N) , * Q » 0 » o 
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3.  THE  METHOD  OF  CHARACTERISTICS 


3(a).  Characteristic  Equations 

Equations  (10)  are  in  general  nonlinear  and  it  is  very  difficult 
to  obtain  closed  form  solutions.  A convenient  method  of  approach  is 
based  upon  obtaining  the  N + 5 base  characteristics  for  this  set  of 
equations;  solution  may  then  be  obtained  by  numerical  integration 
along  these  characteristics.  Proceeding  as  in  [8,9],  we  determine 
the  base  characteristics  as 


(— ) = +v  = + 
W 1 " 


+ 6 + Y - /(l  + 6 + y )2-4y  i 


(11a) 


(f ) = ±v2  = ±1 


N+l 

and  (-^7)  = 0 


(lib) 

(11c) 


Along  the  base  characteristics  dx/dt  = +V2>  we  ^ave 

characteristic  conditions 


(t  " CaVj)(dv-doVj-  l Vjdt)  - f-  deVj-adQV^ 

- - QV2dt  - a2doV^  =0  , j - 1,2  (12a) 

T J J 


and  along  dx/dt  = -V  ,-V2  * we  ^ave  t^ie  conditions 


N 3c . 


(^  -CjV2)(dv+doVj+  l Vjdt)+  ^ dOVj-adQV2-  ^QV2dt+a2davJ=*  0 

, J - 1,2  (12b) 


Also,  along  the  remaining  N+l  base  characteristics  dx/dt  = 0, 
we  require 
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dc 
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, 1 
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(13a) 


(13b) 


In  addition  to  the  above  characteristic  conditions  jump  condi- 
tion must  hold  across  the  wave  fronts,  i.e<,  across  the  base  charac- 
teristics defined  by  dx/dt  = V^,  and  passing  through  the  point 
(0,0).  Using  the  bracket  notation  for  jumps  ([f]=  f -f+)  we  may 
obtain  from  governing  equations  (10) 

N 

I *Ei^j  + = " *^j  (14a) 
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(14b) 

t ] j = [o] j 

(14c) 

[e2^j  = °>  Igfl  = l |o|nsgn(o)] . 

J J J 

(14d) 

“ 0,  t-gplj  - kjtM  sgn(o)]j,  i «=  3 N 

(14e) 

[ eT] j = a 1 6 ] j 

(14f) 

lQ1J-«vJ[a]j-Covj[e]j 
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[Qlj-  ^Wy  iQlj  - -Mflj.  j - 1.2  (14h) 

Note  that  at  the  leading  wave  front  •»  0 since  the  rod  is  initially 
stress  free,  and  thus  the  second  of  equations  (14d,14e)  simplify  to 


3e„ 


l3t~]l  = 


3e 


- k±|  [o]1|  isgnto]1 


(15a) 


i = 3 N (15b) 


From  characteristic  condition  (12a)  and  jump  conditions  (14,15) 
(with  0+  = 0)  , we  obtain  a differential  equation  for  the  stress  jump 
at  the  leading  wave  front  as 
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(16) 


Similarly,  at  the  lagging  wave  front  we  obtain 

tC0v2  y2  dtCTU  N q, 
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3(b).  Boundary  Conditions  and  Initial  Values 

On  the  boundary  x = 0 the  stress  and  temperature  are  prescribed 
by  step  inputs;  i.e.,  o(Q,t)  * H(t)  and  0(0, t)  ■ 6QH(t)  » where 


ft 


12 


H(t)  is  the  unit  step  function.  From  Equation  (10)  we  then  obtain 
on  the  boundary  for  t >.  0 


o(0,t)  - 1 


e^O.t)  - 1 


e2(0,t)  - t , ~<0, t)  - 1 


1 "t/Ti  3ei  "t/Ti 

ei(0,t)  “ "T“(1  " e >•  3t  (0,t)  “ ki  e » 1-3 N 

Pi 1 


eT(0,t)  * a0 
1 o 


0(0, t)  - 0 


Note  that  v and  Q on  the  boundary  cannot  be  found  directly  from 
the  above  prescribed  boundary  condition. 

We  now  seek  to  determine  the  initial  values  of  all  variables  at 
(x,t)  » (0,0)  by  combining  jump  conditions  (14)  with  boundary  condi- 
tions (18) , since  this  point  is  both  on  the  wave  front  and  on  the 
boundary.  We  first  decompose  the  initial  stress  and  temperature  jumps 

into  two  parts  [a  ],,  [0  ],  and  [o  ]„  , [ 0 ] 0 , which  are  assumed 

Ox  ol  O z Ox 

to  propagate  at  the  speeds  and  respectively.  Thus  we  obtain 


l°o]l  + l°o]2  " 1 and  t9o]l  + [6o]2  = 9o 


from  Equations  (18),  and  then  use  Equations  (14g,14h)  to  get 
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Equations  (19,20)  are  then  solved  simultaneously  to  yield  the  follow- 
ing results  at  the  initial  point: 
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(21a) 


(21b) 


(21c) 


(21d) 


From  Equations  (21)  and  the  jump  conditions  (14),  we  obtain  at 
the  leading  wave 


[Vl  “ V1  ^O1! 
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Similarly,  we  obtain  at  the  lagging  wave 


^Vo]2  " " V2  ^°oJ2 
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Finally,  by  combining  Equations  (22)  with  (23)  we  may  then  obtain 
the  desired  initial  values  for  all  variables  at  (x,t)  *(0,0)  as 
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El(0,0)  - 1 


3e„ 


e2(0*0>-°*  bT"1 


3e. 


ei<°«°>  - °>  jrm  ki  * 1 " 3 N 


i 
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i 


eT(0,0) 


6(0,0)  « 6 


Q(0,0)  - + 1^2) 

T Vx  V2 


(24) 


With  characteristic  conditions  (12 ,13) , jump  conditions  (14), 
boundary  condition  (18)  and  initial  values  (24) , a numerical  procedure 
for  obtaining  stress  and  temperature  solutions  may  be  established. 
However,  closed  form  solutions  may  be  obtained  in  some  special  cases, 
and  we  first  develop  such  a solution  in  the  next  section  for  the 
uncoupled  thermal  wave  problem.  We  also  solve  this  special  case  by 
the  method  of  characteristics,  before  proceeding  to  the  numerical 
solution  of  the  general  problem. 


I 


» 


' 


i 

> 


rT~  wr»m 


is 


4.  UNCOUPLED  THERMAL  WAVE  PROBLEM 

4(a).  Closed  Form  Solution 

For  the  uncoupled  case  we  set  a = 0 in  Equation  (lOf)  and  obtain 
the  uncoupled  energy  equation 


3x 


(25) 


By  combining  Equations  (10a)  with  (25)  we  may  then  obtain  an  uncoupled 
thermal  wave  equation  as 


a2e  l _ae 

at2  T 3t 


2 

.2  3 0 
T 

3x 


0 


(26) 


where  V_  ■ /k/C  t is  the  uncoupled  thermal  wave  speed.  (Note  that 
T o 

VT  = Lim  V9  via  Equation  (6)  while  the  uncoupled  mechanical  wave 
T cr*0 

speed  V ■=  Lim  V1 . ) We  can  simplify  the  problem  further  by  intro- 
6 a-X) 

ducing  a new  set  of  coordinates  as 


(27) 


Equation  (26)  then  becomes 


a2e  ae_  a2e 

*2  + * *2 

at  at  ax 


(28) 


The  rod  is  initially  at  uniform  temperature  Tq  , and  a step 

* * + 

in  temperature  is  applied  on  the  boundary  x * 0 at  t * 0 
Thus,  the  initial  values  are  given  in  this  case  as 


1 


/ 


■ 
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6(x  , 0 ) - 0 


'1 


— *(x*.  0")  - 0 
3t 


and  the  boundary  condition  as 


6(0, t*)  - 6 H(t*) 


(29a) 


(29b) 


With  the  use  of  the  Laplace  transform  of  Equations  (28,29),  we 
may  obtain  the  solution  for  the  temperature  distribution  along  the 
rod  as 


e«e 


.{f 


* x* 


i t'^I1  (1/t  2 -x*  2 ) H ( t ' -x*  ) + j e-t,^H(t'-x*) 


dt' 


+ e 


Al(t  *-**>} 


(30) 


where  I^(t)  is  the  modified  Bessel  function  of  the  first  kind  of 

it  it 

order  one.  At  the  wave  front,  i.e.,  at  x ■ t , Equation  (30) 


reduces  to 
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which  can  be  readily  evaluated  by  numerical  integration. 

4(b).  Numerical  Solution  via  Method  of  Characteristics 

Setting  u - 0 in  Equations  (12)  we  obtain  the  uncoupled  mech- 
anical characteristic  conditions  (see  [5]) 


N 3e 

dv  *=  ±(do  + £ -r—  dt) 

i-2  n 


along  dx/dt  * ±Vfi  respectively , and  the  uncoupled  thermal  charac- 
teristic conditions 


d0  - +(tdQ  + Qdt)^—  (34) 

along  dx/dt  - ±VT  respectively.  The  solution  to  the  uncoupled 
mechanical  wave  problem  by  the  method  of  characteristics  has  been  dis- 
cussed in  detail  by  Cozzarelli  and  Shaw  [8-10] , and  thus  we  will 
confine  our  attention  to  the  numerical  solution  of  the  uncoupled 
thermal  wave  problem. 

Equations  (34)  are  firs4-  as 


d0  - +(dQ*  + Q*  dt*) 


along  dx  /dt  ■ ±1  respectively,  where 


* QV 


Next  we  note  that  jump  condition  (14h)  simplifies  to 
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* 

Q - 6 


Combining  Equations  (35,37)  we  obtain  at  the  uncoupled  thermal  wave 
front 


2d0  - -0dt 


This  result  Integrates  Immediately  to 


0 - 0 e 
o 


which  is  Identical  with  Equation  (31).  Next,  the  boundary  condition 
* 

for  t 0 Is 


0(O,t  ) « 0 


Finally,  the  Initial  values  are  given  by 


0(0,0)  = 0, 


Q (0,0)  - 0 


The  equations  developed  in  this  section  may  be  approximated  by 

a finite  difference  scheme.  We  use  mesh  points  at  the  intersections 

* * 

of  the  three  base  characteristic  families  dx  /dt  ■ +1,  -1  and  0 , 

it 

which  gives  a constant  increment  Ax  ■ At  . This  is  shown  in 
Figure  1.  Mesh  points  are  distinguished  as  wave  front  points,  i.e., 

* it 

lying  along  the  wave  front  (x  ■ t ),  boundary  points,  i.e.. 


I 


lying  on  the  boundary  (x  = 0) , and  interior  points 


We  use  SAt  to  indicate  the  discrete  values  of  t at  which 


calculation  are  performed,  and  (J-l)Ax  to  indicate  the  discrete 


1 representing  the  boundary  X **  0»  Because 


values  of  X with  J 


of  the  alternating  pattern  of  mesh  points  defined  by  the  intersections 


of  the  three  base  characteristics,  one  row  of  constant  t mesh 


points  will  start  from  the  boundary  (J  = 1)  with  an  increment  of 


2Ax  in  the  x direction  until  the  wave  front  is  reached,  using 


only  odd  values  of  J . The  next,  row  will  start  at  a taesh  point 


(J  *■  2)  with  ah  increment  of  2 Ax  until  the  wave  ftont  is  reached 


using  only  even  values  of  J . It  is  convenient  from  a computer 


storage  point  of  view  to  condense  these  two  rows  such  that  one  value 


will  contain  all  values 


of  K , which  determines  two  values  of  t 


This  is  done  by  letting 


only  odd  Values  of  J 


(2K+l)At  only  even  values  of  J 


Ordinary  differential  equations  (35)  may  be  approximated  by  finite 
difference  equations  using  a modified  Euler  method  (see  [11]),  thus 


The  resulting  algebraic  equations  at  the  interior  points  are  given  by 
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0 (J,S)-0 (J-l , S-l)=  -Q*(J,S)+  Q*(J-1,S-1)-  y[Q*(J,S)+  Q*(J-1,S-1) ]At* 


along  dx  /dt  ■ +1,  and 


(44a) 


0(J,S)-0(J-1  » S-l)»  Q*(J,S)-  Q*(J-1,S-1)+  j[Q*(J,S)+  Q* (J-l ,S-1) ] At* 

(44b) 

* * 

along  dx  /dt  = -1. 

Equations  (44)  may  be  solved  for  0(J,S)  and  Q (J,S)  once  the  pre- 
vious (S-l)  values  are  known. 

For  points  on  the  boundary,  J ■ 1 , the  boundary  condition  (40) 

will  be  used  instead  of  Equation  (44a).  Similarly,  for  points  on  the 

wave  front  Equation  (44b)  is  replaced  by  jump  condition  (37).  In 

any  case,  the  initial  values  (41)  must  be  applied  to  start  the  numeri- 

* * 

cal  calculation.  The  temperature  difference  field  6(x  ,t  ) was 

calculated  both  by  the  numerical  procedure  outlined  in  this  section 

and  by  the  closed  form  solution  of  the  previous  section  (Equation 

(30)).  Identical  results  were  obtained  via  the  two  methods,  and 

Figure  2 shows  some  field  solutioiB  for  0 plotted  as  a function  of 

* 

distance  at  several  specific  times  (t  = 1,2,3).  This  figure  clearly 
shows  that  the  use  of  the  modified  Fourier  heat  conduction  law  results 
in  a damped  thermal  wave,  even  when  damping  due  to  viscoelastic  and 
coupling  effects  are  not  present. 


at 
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5.  NUMERICAL  PROCEDURE  AND  RESULTS 
FOR  THE  GENERAL  PROBLEM 

In  the  general  coupled  case  (a  4 0),  the  coupling  of  the  charac- 
teristic condition  (12)  necessitates  a more  complicated  numerical  pro- 
cedure than  in  section  4(b).  However,  the  method  of  approach  for 
a 4 0 contains  some  similarities  to  that  for  the  uncoupled  case. 

5(a).  Numerical  Procedure 

We  choose  - 3 as  an  illustrative  example,  and  use  a finite 
difference  scheme  with  mesh  points  lying  on  the  intersections  of  the 
base  characteristics  of  the  three  families  dx/dt  * +3,  -3  and  0 . 
This  is  shown  in  Figure  3.  Additional  points  along  the  base  charac- 
teristics dx/dt  « +1  and  -1  may  be  obtained  by  linear  interpolation 
between  the  mesh  points,  since  these  additional  points  lie  within  the 
domain  of  dependence  of  the  differential  equatioi.s.  Mesh  points  in 
Figure  3 are  distinguished  as  leading  wave  front  points,  i.e.,  lying 
along  the  leading  wave  front  (x  = 3t),  boundary  points,  i.e.,  lying 
on  the  boundary  (x  « 0),  lagging  wave  front  points,  i.e.,  lying  along 
the  lagging  wave  front  (x  = t)  , and  interior  points  (two  types). 
These  different  types  of  mesh  points  require  separate  calculation 
routines.  We  use  SAt  to  indicate  the  discrete  values  of  t at 
which  calculation  are  performed,  and  (J-l)Ax  to  indicate  the  dis- 
crete values  of  x with  J * 1 representing  the  boundary  x “ 0. 

A similar  numerical  approach  but  with  a cruder  finite  difference  mesh 
scheme  was  used  by  Lopez  and  Lord  [12]  to  study  the  special  case  of 
coupled  linear  thermoelastic  wave  propagation  with  second  sound. 


23 


As  in  the  uncoupled  thermal  wave  problem  we  may  take  advantage 
of  the  alternating  pattern  of  the  finite  difference  mesh  to  save 
computer  storage  space.  We  note  that  one  row  of  constant  t mesh 
points  starts  from  the  boundary  (J=l)  with  an  increment  of  6Ax  in  the 
x direction  up  to  the  leading  wave  front,  using  only  odd  values  of. 

J . The  next  row  starts  at  a mesh  point  (J  = 4)  with  an  increment 
of  6Ax  up  to  the  leading  wave  front,  using  only  even  values  of  J . 

We  may  condense  these  two  rows  by  defining  S in  terms  of  K (see 
Figure  3)  as 


t = SAt 


f2KAt  only  odd  values  of  J 


(2K+l)At  only  even  values  of  J , K = 0,1,2,...  (45) 


Ordinary  differential  equations  (12,13)  may  be  approximated  by 
finite  difference  equations  using  a modified  Euler  method,  i.e. , we 
employ  average  values  for  the  functions  Se^/St  , 3e^,/3t  and  Q . 
The  resulting  nonlinear  algebraic  equations  obtained  from  equations 
(12)  for  interior  points  are  given  by 


<7  - 9C  ) 
t o 


r 3e 


-3At 


v(J ,S)-v(J-3 ,S-l)-3o (J ,S)+3o (J-3 ,S-1) 


2 9e2  9ei  9ei 

7~(J,S)+  -jj-j^-(J-3,S-l)  N -^(J,S)+  — (J-3.S-1) 

_ + £ « 

L i-3  ■ 


- — [0 (J,S)-0 (J-3 ,S-1) ]-9a[Q(J ,S)-Q(J-3,S-1) ] 

T 


At  [Q.(J>s)+9(J-3is-U].  27a2[a(J,S)-a(J-3,S-l)]»  0 


(46a) 


— V- 
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along  dx/dt  = +3,  and  similarly 


W 

(-  - 9Cq)-  v(J,S)-v(J+3,S-1)+3o(J,S)-30(J+3,S-1) 

Sty  9fc2  dt ' 

— -(J,S)+  — (J+3.S-1)  N ~(J,S)+  —(J+3.S-1) 

+3At  

L 7 i*=3  2 


+—  [0(J,S)-0(J+3,S-1)]  - 9a[Q(J,S)-Q(J+3,S-l)] 


- 5a  At 

T 


-]+  27a  [a(J,S)-o(J+3,S-l>]  = 0 


(46b) 


along  dx/dt  = -3  . Also, 


- Cq) jv(J,S)-v(J-l,S-l)-o(J,S)+a(J-l,S-l) 

- 3e2  3eo  ■ 3e 

— (J,s>+  — ^(J-1,S-D  N 7p(J,S)+r^(J-l,S-l) 

-At  ^ f + l 

L i=3 


- [0 (J,S)-0 (J-1,S-1) ]-a[Q(J ,S)-Q(J-1,S-1) ] 


a A r ( 
-TAt[- 


-]  -a  [a(J,S)-o(J-l,S-l)]  = 0 


(46c) 


along  dx/dt  = +1,  and 


(t  "Ca){V ( J ’ S) ~v ( J+1 » S_1 )+o (J , S) -a ( J+l , S-l ) 

' de2  3e?  3e-»  3e, 

T— (J,S)+ rr£(J+1,s-1)  N — -(J,s)+  — i(j+i,s-i) 

+At  — — + y 

2 L 7 

L i=3  1 


+ ~ [8 (J ,S)-0 (J+l , S-l) ]-a[Q(J ,S)-Q(J+1 ,S-1) ] 
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- At  >g.)tQ(J+l  >S~Uj+  a2[o(J,S)-o(J+l,S-l)]-  0 
T i 


(46d) 


along  dx/dt  = -1.  Finally,  equations  (13)  yield 


e1(J,S)-e1(J,S-2) 


e2(J,S)-E2(J,S-2) 


e.(J,S)-e1(J,S-2) 


eT(J,S)-eT(J,S-2) 


2 St  ?r<J-s)+  tt<j-s-2) 


"3  € ry  3 £ « 


n>q  5ei  1 

2 At  ^r(J’S)+  ir(J»S-2)  , i=3 , . 


*eT  3eT 

— (J.SH  ■g~(J«S-2) 


(46e) 


(46f) 


. ,N  (46g) 


(46h) 


along  dx/dt  = 0 . Nonlinear  algebraic  equations  (46a-h)  may  be  solved 
for  v,  o,  e^,  Se^/St,  e,  Q,  e , and  3eT/3t  for  mesh  points  at  a 
given  S once  the  previous  (S-l)  values  are  known. 

The  methods  used  to  obtain  solutions  at  the  various  types  of 
mesh  points  shown  in  Figure  3 are  summarized  below. 

(i)  Initial  Point:  All  initial  values  are  given  by  Equations 
(21-24). 

(ii)  Leading  Wave  Front  Points:  Stress  is  obtained  by  integ- 
rating Equation  (16),  and  the  remaining  variables  are  obtained  from 
Equations  (14). 


(iii)First  Boundary  Point  (see  Figure  4a):  At  the  initial  point 
(0,0)  all  N + 5 variables  are  known  from  Equations  (24).  At  the 
first  boundary  point  (0,2)  we  know  the  values  of  all  variables 
except  v and  Q from  Equations  (18).  We  start  an  iteration  process 
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* 


1 


for  determining  the  values  of  v and  Q by  first  setting  v and  1} 
at  point  (0,2)  equal  to  their  values  at  point  (0,0).  We  interpolate 
for  values  of  any  variable  at  point  1 using  the  values  at  points 

(0,0)  and  (0,2)  . At  point  4 all  variables  are  known  from  the  previous 

time  step.  The  jumps  of  all  variables  at  points  2 and  6 are  obtained 
from  the  integration  of  Equation  (17)  and  jump  conditions  (14)  (see 
enlarged  Figure  4b).  Basing  on  linear  interpolation  and  the  assumption 
that  the  slope  of  any  variable  at  the  lagging  wave  front  will  remain 
the  same  as  it  crosses  the  wave  front,  we  obtain 

©=(®  - [2]  - ®)  x | + © + 12] 

(D=((D  “ 12]  - ©)  x 5 + © 

©-(©  - 12]  - ©)  x \ + © (47) 

where  ® indicates  the  value  of  a variable  at  point  1 , and 
[2]  represents  the  jump  across  the  lagging  wave  front  at  point  2. 

Once  all  values  of  the  variables  at  points  (1-5)  are  obtained,  we  may 
use  Equations  (14,46)  to  solve  for  the  N + 5 variables  at  point  6. 

From  points  2 and  6 and  Equations  (46b, 46d)  we  may  then  find  an 
improved  value  of  v and  Q at  point  (0,2).  By  this  procedure  we 
may  obtain  v and  Q at  point  (0,2)  through  iteration 

(iv)  Subsequent  Boundary  Points  (see  Figure  5):  First  we  assume 
values  for  v and  Q at  point  (0,S).  Then  we  interpolate  point  1 

from  points  (0,S-2)  and  (0,S).  We  know  point  3 from  previous  time 
steps.  Therefore,  point  2 can  be  obtained  from  points  1 and  3 via 
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linear  interpolation.  Finally  we  may  evaluate  the  values  of  v and 
Q at  point  (0,S)  from  points  2 and  3 from  Equations  (46b, 46d)  and 
iteration  as  in  part  (iii). 

(v)  Lagging  Wave  Front  Points  (see  Figure  6) : The  values  of 
the  N + 5 variables  at  points  1,  4 and  5 are  known  from  previous 
time  steps.  The  jumps  at  points  2 and  6 can  be  obtained  from  Equations 
(14,17).  The  values  of  the  variables  at  points  2 and  3 can  be  found 

by  linear  interpolation.  The  information  needed  to  evaluate  point  6 
from  equations  (46)  is  then  complete. 

(vi)  Interior  Points  - Type  A (see  Figure  7):  We  know  the 
values  of  all  variables  at  points  1,  7 and  9 from  previous  time  steps, 
and  the  jumps  at  points  4 and  8 are  obtained  from  Equations  (14,17). 

We  then  use  linear  interpolation  to  obtain  points  2,  3,  4,  5 and  6. 

Point  8 may  then  be  obtained  from  points  3,  4,  5,  6 and  7 as  in 

part  (v).  Once  points  1,  2,  4,  8 and  are  known,  we  may  evaluate 
point  10  from  Equations  (46). 

(vii)  Interior  Points  - Type  B (see  Figure  8) : The  values  of  all 
variables  are  known  at  points  1,  4 and  5,  and  points  2,  3 can  be  obtained 
from  linear  interpolation.  Once  points  1,  2,  3,  4 and  5 are  known, 
we  then  have  all  information  necessary  to  solve  for  point  6. 

In  the  next  section  we  present  some  numerical  results  with 
= 3 for  several  illustrative  cases,  as  obtained  from  a computer 
program  written  for  the  CYBER  173.  For  arbitrary  the  lagging  wave 

front  point  is  not  necessarily  a mesh  point,  but  values  at  the  point 
may  still  be  easily  calculated  via  linear  interpolation. 


5(b).  Discussion  of  Results 

Although  the  program  written  can  be  employed  for  various  gener- 
alized Kelvin  material  models  with  up  to  twelve  elements  (i.e.,  one 
linear  elastic  spring,  one  nonlinear  dashpot  and  ten  nonlinear  Kelvin 
elements),  only  three  material  models  will  be  discussed  here.  These 
models  are  the  linear  Maxwell  model  (n  = 1) , a nonlinear  Maxwell 
model  (n  = 3)  and  the  linear  Maxwell-Kelvin  model.  These  three  models 
are  sufficient  to  illustrate  the  effects  of  second  sound,  thermo- 
mechanical coupling  and  nonlinear  viscoelastic  damping.  Also,  all 
models  will  be  solved  with  a positive  applied  step  stress  at  the  boun- 
dary, but  with  a zero  temperature  increment  (i.e.,  0o=O) . An  applied 
constant  tensile  stress  will  cause  the  temperature  to  tend  to  drop  at 
the  boundary,  and  thus  heat  must  be  continuously  added  at  x = 0 to 
maintain  6 at  zero. 

In  selecting  a set  of  illustrative  values  for  the  nondimens ional 

material  constants  k,  t,  C and  a , we  seek  a case  in  which  coupling 

a 

and  second  sound  are  significant  effects.  However,  on  physical  grounds 
we  would  still  expect  in  general  that  in  heat  conduction  equation  (10a) 
k and  also  that  in  energy  equation  (lOf)  Ca-|~  > a . 

dX  at  dt  dt 

Thus  we  will  choose  k > t and  > a , where  this  latter  condition 
also  ensures  non-negative  values  for  6 and  y [see  Equation  (9b)]. 
Furthermore,  we  must  also  satisfy  conditions  (6b)  which  in  nondimen- 
sional  form  becomes  for  the  present  case  = 3 > = 1 , 

thereby  ensuring  that  the  leading  wave  be  predominantly  thermal.  In 
consideration  of  the  above  we  select  k = 1.0,  t k 0.3  , = 0.7 

and  a = 0.5  , which  then  yields  6 ■ 0.555,  y * 5.818,  VT  * 2.182 
and  Vg  = 1.273  via  Equations  (9b,  11a).  All  numerical  results  obtained 
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using  these  values  were  checked  by  substitution  into  the  governing 
equations  with  the  use  of  the  five-point  finite  difference  derivative 


formula 


{o  ■ 1 k (f-2-  8,-l  + 8fl  - f2>  + 35  £V<{) 


where  f^  is  the  first  derivative  of  f at  the  central  point,  h 
is  the  increment,  and  the  last  term  is  the  truncation  error. 

For  the  first  case  of  a linear  Maxwell  model.  Figures  9 and  10 
give  the  field  solutions  for  the  nondimens ional  stress  and  temperature 
difference  respectively  as  functions  of  distance  at  times  t * 0.5,  1.5. 
We  see  from  these  two  figures  that  second  sound  results  in  two  wave 
fronts  with  finite  speeds  along  which  both  the  mechanical  and  thermal 
disturbances  propagate.  Note  that  the  stress  (Figure  9)  decreases 
monotonically  from  the  boundary  value  of  1.0  to  a value  behind  the  lag- 
ging wave  front.  After  a negative  jump  the  stress  again  decreases 
monotonically  to  a value  behind  the  leading  wave  front.  Upon  compar- 
ing these  stresses  with  similar  results  from  uncoupled  theory  in  P-0] , 
we  find  that  thermomechanical  coupling  is  a damping  effect  which 
results  in  a faster  relaxation  of  the  stress  jumps.  The  temperature 
difference  (Figure  10)  decreases  from  the  prescribed  boundary  value 
to  a value  behind  the  lagging  wave  front.  After  a negative  jump  at 
the  lagging  wave  front  it  increases  to  a value  at  the  leading  wave 
front,  which  is  less  than  the  initial  value  of  zero.  The  temperature 
difference  solution  in  the  uncoupled  case  for  this  problem  is  zero 
everywhere. 

These  figures  show  that  the  damping  effect  due  to  thermomechanical 
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of  all  jumps  at  the  two  wave  fronts.  The  stress  jump  at  the  leading 
wave  front  is  much  smaller  than  that  at  the  lagging  wave  front,  which 
is  as  expected  for  a mechanical  input.  The  difference  in  magnitude 
between  the  two  temperature  jumps  is  not  significant,  but  we  see  that 
they  are  of  opposite  sign.  Also  we  see  that  at  the  lagging  wave  front 
the  stress  and  temperature  jumps  are  of  the  same  sign,  while  at  the 
leading  wave  front  they  are  of  opposite  sign.  In  general.  Figure  10 
indicates  that  the  bulk  of  the  thermal  energy  is  carried  by  the  lagging 
wave  front  for  a mechanical  input  of  this  type.  A similar  result  was 
obtained  in  [12]  for  a linear  thermoelastic  material. 

Figures  11  and  12  give  the  results  for  the  nonlinear  Maxwell 
model  with  n = 3.  On  comparing  Figures  9 and  11  we  see  that  increas- 
ing n results  in  a slower  relaxation  of  the  stress  field.  This 
follows  from  constitutive  equation  (10c),  since  as  n increases,  with 
| a ] il,  the  effective  viscosity  of  the  dashpot  increases.  More  energy 
will  be  expended  in  displacing  a dashpot  as  n increases,  and  accord- 
ingly Figure  12  shows  a greater  temperature  drop  than  does  Figure  10. 

Finally,  the  solutions  for  the  linear  Maxwell-Kelvin  model  with 
x ■ y = q - 1 are  shown  in  Figures  13  and  14.  We  note  by  com- 
parison with  Figure  9 that  adding  a linear  Kelvin  element  to  a linear 
Maxwell  element  results  in  a faster  relaxation  of  the  stress  field. 

This  is  clearly  due  to  the  additional  viscoelastic  damping  caused  by 
the  Kelvin  element.  This  damping  effect  also  results  in  a smaller 
temperature  drop  than  for  the  linear  Maxwell  model,  as  Bhovn  in  Figure 
14. 
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6 . BRIEF  SUMMARY 

In  order  to  ensure  a finite  speed  of  propagation  for  a thermal 
disturbance,  a modified  Fourier  heat  conduction  law  was  used.  A one- 
dimensional single  integral  nonlinear  thermoviscoelastic  constitutive 
relation  and  a linearized  thermomechanically  coupled  energy  equation 
were  also  employed.  Using  these  relations  together  with  the  one- 
dimensional strain-displacement  relation  and  the  equation  of  motion, 
the  one-dimensional  wave  propagation  problem  was  then  studied  for  a 
semi-infinite  rod. 

Due  to  the  nonlinearity  and  general  complexity  of  the  governing 
equations,  a numerical  approximation  based  on  the  method  of  character- 
istics was  employed  to  obtain  the  field  solutions  for  stress  and  tem- 
perature. Two  thermomechanically  coupled  wave  fronts  with  finite  speeds 
were  generated  along  the  positive  axis  by  applied  step  stress  and 
temperature  inputs  at  the  boundary.  A closed-form  solution  in  integral 
form  was  obtained  via  the  method  of  Laplace  Transformation  for  the 
special  case  of  an  uncoupled  thermal  wave.  A numerical  solution  based 
on  the  method  of  characteristics  was  also  obtained  in  this  case. 

For  the  general  coupled  case,  the  initial  values  were  obtained 
from  a decomposition  of  the  initial  input  and  from  the  jump  conditions. 

A finite  difference  approximation  using  the  modified  Euler  method  was 
used  for  the  nonlinear  characteristic  equations.  Then  the  field  solu- 
tions were  obtained  from  these  nonlinear  algebraic  equations  via 
Iteration.  A detailed  discussion  of  the  numerical  procedure  was  given. 
Numerical  results  were  given  and  discussed  for  three  viscoelastic 


models. 
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